Stress x EE MA
Loading packages
knitr::opts_chunk$set(echo = TRUE)
pacman::p_load(tidyverse,
here,
metafor,
emmeans,
orchaRd,
MCMCglmm,
corpcor,
clubSandwich,
MuMIn,
patchwork,
naniar,
GoodmanKruskal,
networkD3,
ggplot2,
ggalluvial,
ggthemr,
cowplot)
# needed for model selection using MuMin within metafor
eval(metafor:::.MuMIn)Loading data and functions
dat <- read_csv(here("Data","Full_data.csv"))
# Load custom function to extract data
#source(here("R/functions.R"))
source(here("R/functions_cleanedup.R")) Data organisation
Removing study with negative values, getting effect sizes from function, ‘flipping’ effect sizes so that all effect sizes are higher values = individuals do better and learning/memory, shifting negative values to positive as lnRR cannot use negative values, assigining human readable terms, and creating VCV of variance
#removing study with negative values as these are unable to be used for lnRR
dat <- droplevels(dat[!dat$First_author == 'Wang',])
#Getting effect sizes
effect_size <- effect_set(CC_n = "CC_n", CC_mean = "CC_mean", CC_SD = "CC_SD",
EC_n = "EC_n", EC_mean = "EC_mean" , EC_SD ="EC_SD",
CS_n = "CS_n", CS_mean = "CS_mean", CS_SD = "CS_SD",
ES_n = "ES_n", ES_mean = "ES_mean", ES_SD = "ES_SD",
data = dat)
#'pure' effect sizes
effect_size2 <- effect_set2(CC_n = "CC_n", CC_mean = "CC_mean", CC_SD = "CC_SD",
EC_n = "EC_n", EC_mean = "EC_mean" , EC_SD ="EC_SD",
CS_n = "CS_n", CS_mean = "CS_mean", CS_SD = "CS_SD",
ES_n = "ES_n", ES_mean = "ES_mean", ES_SD = "ES_SD",
data = dat)
#Removing missing effect sizes
dim(dat)
full_info <- which(complete.cases(effect_size) == TRUE)
dat_effect <- cbind(dat, effect_size, effect_size2)
dat <- dat_effect[full_info, ]
names(dat_effect)
dat <- dat_effect[full_info, ]
dim(dat_effect)
dimentions <- dim(dat)
#Flipping 'lower is better' effect sizes
#flipping lnRR for values where higher = worse
dat$lnRR_Ea <- ifelse(dat$Response_direction == 2, dat$lnRR_E*-1,ifelse(is.na(dat$Response_direction) == TRUE, NA, dat$lnRR_E)) # currently NAswhich causes error
dat$lnRR_Sa <- ifelse(dat$Response_direction == 2, dat$lnRR_S*-1,ifelse(is.na(dat$Response_direction) == TRUE, NA, dat$lnRR_S)) # currently NAswhich causes error
dat$lnRR_ESa <- ifelse(dat$Response_direction == 2, dat$lnRR_ES*-1,ifelse(is.na(dat$Response_direction) == TRUE, NA, dat$lnRR_ES)) # currently NAswhich causes error
#flipping 'pure effect sizes'
dat$lnRR_E2a <- ifelse(dat$Response_direction == 2, dat$lnRR_E2*-1,ifelse(is.na(dat$Response_direction) == TRUE, NA, dat$lnRR_E2)) # currently NAswhich causes error
dat$lnRR_S2a <- ifelse(dat$Response_direction == 2, dat$lnRR_S2*-1,ifelse(is.na(dat$Response_direction) == TRUE, NA, dat$lnRR_S2)) # currently NAswhich causes error
dat$lnRR_ES2a <- ifelse(dat$Response_direction == 2, dat$lnRR_ES2*-1,ifelse(is.na(dat$Response_direction) == TRUE, NA, dat$lnRR_ES2)) # currently NAswhich causes error
dat$lnRR_E3a <- ifelse(dat$Response_direction == 2, dat$lnRR_E3*-1,ifelse(is.na(dat$Response_direction) == TRUE, NA, dat$lnRR_E3)) # currently NAswhich causes error
dat$lnRR_S3a <- ifelse(dat$Response_direction == 2, dat$lnRR_S3*-1,ifelse(is.na(dat$Response_direction) == TRUE, NA, dat$lnRR_S3)) # currently NAswhich causes error
#flipping SMD
dat$SMD_Ea <- ifelse(dat$Response_direction == 2, dat$SMD_E*-1,ifelse(is.na(dat$Response_direction) == TRUE, NA, dat$SMD_E)) # currently NAswhich causes error
dat$SMD_Sa <- ifelse(dat$Response_direction == 2, dat$SMD_S*-1,ifelse(is.na(dat$Response_direction) == TRUE, NA, dat$SMD_S)) # currently NAswhich causes error
dat$SMD_ESa <- ifelse(dat$Response_direction == 2, dat$SMD_ES*-1,ifelse(is.na(dat$Response_direction) == TRUE, NA, dat$SMD_ES))
#rounding down integer sample sizes
dat$CC_n <- floor(dat$CC_n)
dat$EC_n <- floor(dat$EC_n)
dat$CS_n <- floor(dat$CS_n)
dat$ES_n <- floor(dat$CS_n)
#assigning human readable terms
dat <- dat %>% mutate(Type_learning = case_when(Type_learning == 1 ~ "Habituation",
Type_learning == 2 ~ "Conditioning",
Type_learning == 3 ~ "Recognition",
Type_learning == 4 ~ "Unclear"),
Learning_vs_memory = case_when(Learning_vs_memory == 1 ~ "Learning",
Learning_vs_memory == 2 ~ "Memory",
Learning_vs_memory == 1 ~ "Unclear"),
Type_reinforcement = case_when(Type_reinforcement== 1 ~"Appetitive",
Type_reinforcement== 2 ~ "Aversive",
Type_reinforcement== 3 ~ "Not applicable",
Type_reinforcement== 4 ~ "Unclear"),
Type_stress_exposure = case_when(Type_stress_exposure == 1 ~ "Density",
Type_stress_exposure == 2 ~ "Scent",
Type_stress_exposure == 3 ~ "Shock",
Type_stress_exposure == 4 ~ "Exertion",
Type_stress_exposure == 5 ~ "Restraint",
Type_stress_exposure == 6 ~ "MS",
Type_stress_exposure == 7 ~ "Circadian rhythm",
Type_stress_exposure == 8 ~ "Noise",
Type_stress_exposure == 9 ~ "Other",
Type_stress_exposure == 10 ~ "Combination",
Type_stress_exposure == 11 ~ "unclear"),
Age_stress_exposure = case_when(Age_stress_exposure == 1 ~ "Prenatal",
Age_stress_exposure == 2 ~ "Early postnatal",
Age_stress_exposure == 3 ~ "Adolescent",
Age_stress_exposure == 4 ~ "Adult",
Age_stress_exposure == 5 ~ "Unclear"),
Stress_duration = case_when(Stress_duration == 1 ~ "Acute",
Stress_duration == 2 ~ "Chronic",
Stress_duration == 3 ~ "Intermittent",
Stress_duration == 4 ~ "Unclear"),
EE_social = case_when(EE_social == 1 ~ "Social",
EE_social== 2 ~ "Non-social",
EE_social == 3 ~ "Unclear"),
EE_exercise = case_when(EE_exercise == 1 ~ "Exercise",
EE_exercise == 2 ~ "No exercise"),
Age_EE_exposure = case_when(Age_EE_exposure == 1 ~ "Prenatal",
Age_EE_exposure == 2 ~ "Early postnatal",
Age_EE_exposure == 3 ~ "Adolescent",
Age_EE_exposure == 4 ~ "Adult",
Age_EE_exposure == 5 ~ "Unclear"),
Exposure_order = case_when(Exposure_order == 1 ~ "Stress first",
Exposure_order == 2 ~ "Enrichment first",
Exposure_order == 3 ~ "Concurrently",
Exposure_order == 4 ~ "Unclear"),
Age_assay = case_when(Age_assay == 1 ~ "Early postnatal",
Age_assay == 2 ~ "Adolescent",
Age_assay == 3 ~ "Adult",
Age_assay == 4 ~ "Unclear"),
Sex = case_when(Sex == 1 ~ "Female",
Sex == 2 ~ "Male",
Sex == 3 ~ "Mixed",
Sex == 4 ~ "Unclear"))
#making variance VCVs
VCV_E <- impute_covariance_matrix(vi = dat$lnRRV_E, cluster = dat$Study_ID, r = 0.5)
VCV_S <- impute_covariance_matrix(vi = dat$lnRRV_S, cluster = dat$Study_ID, r = 0.5)
VCV_ES <- impute_covariance_matrix(vi = dat$lnRRV_ES, cluster = dat$Study_ID, r = 0.5)
VCV_Ea <- impute_covariance_matrix(vi = dat$SMDV_E, cluster = dat$Study_ID, r = 0.5)
VCV_Sa <- impute_covariance_matrix(vi = dat$SMDV_S, cluster = dat$Study_ID, r = 0.5)
VCV_ESa <- impute_covariance_matrix(vi = dat$SMDV_ES, cluster = dat$Study_ID, r = 0.5)Data exploration
General
#Number of effect sizes
length(unique(dat$ES_ID)) ## [1] 92#Number of studies
length(unique(dat$Study_ID))## [1] 30#Publication years
min(dat$Year_published) ## [1] 2006max(dat$Year_published)## [1] 2021Explore associations among predictor variables
#see columns with missing values
vis_miss(dat) +
theme(plot.title = element_text(hjust = 0.5, vjust = 3),
plot.margin = margin(t = 0.5, r = 3, b = 1, l = 1, unit = "cm")) +
ggtitle("Missing data in the selected predictors") #no mising values
#useGoodman and Kruskal’s τ measure of association between categorical predictor variables (function from package GoodmanKruskal: https://cran.r-project.org/web/packages/GoodmanKruskal/vignettes/GoodmanKruskal.html)
#GKmatrix <- GKtauDataframe(subset(dat, select = c("Sex", "Type_learning", "Learning_vs_memory", #"Type_reinforcement", "Type_stress_exposure", "Age_stress_exposure", "Stress_duration", #"EE_social_HR", "EE_exercise", "Age_EE_exposure", "Exposure_order", "Age_assay")))
#plot(GKmatrix)
#simple pairwise contingency tables
table(dat$Type_learning, dat$Type_reinforcement)
table(dat$Age_stress_exposure, dat$Age_EE_exposure)
table(dat$Type_stress_exposure, dat$Age_stress_exposure)
table(dat$Type_stress_exposure, dat$Age_assay)
table(dat$Type_stress_exposure, dat$Stress_duration)Alluvial diagrams
Species and strain
# create a frequency table for Species and Strain variables
freq_1 <- as.data.frame(table(dat$Common_species, dat$Strain)) %>% rename(Species = Var1, Strain = Var2)
is_alluvia_form(as.data.frame(freq_1), axes = 1:2, silent = TRUE)
ggplot(data = freq_1,
aes(axis1 = Species, axis2 = Strain, y = Freq)) +
geom_alluvium(aes(fill = Strain)) +
geom_stratum(aes(fill = Species)) +
geom_text(stat = "stratum", aes(label = after_stat(stratum))) +
#theme_minimal() +
theme_void() +
theme(legend.position = "none",
plot.title = element_text(hjust = 0.5, vjust = 3),
axis.title.x = element_text(),
axis.text.x = element_text(face="bold")) +
scale_x_discrete(limits = c("Species", "Strain"), expand = c(0.15, 0.05), position = "top") +
ggtitle("Species and strains used (set1)")
Species, strain and sex
# create frequency table for Sex, Species and Strain
freq_2 <- as.data.frame(table(dat$Sex, dat$Common_species, dat$Strain)) %>% rename(Sex = Var1, Species = Var2, Strain = Var3)
is_alluvia_form(as.data.frame(freq_2), axes = 1:3, silent = TRUE)
ggplot(data = freq_2,
aes(axis1 = Sex, axis2 = Species, axis3 = Strain, y = Freq)) +
geom_alluvium(aes(fill = Sex)) +
geom_flow() +
geom_stratum(aes(fill = Sex)) +
geom_text(stat = "stratum", aes(label = after_stat(stratum))) +
#theme_minimal() +
theme_void() +
theme(legend.position = "none",
plot.title = element_text(hjust = 0.5, vjust = 3),
axis.title.x = element_text(),
axis.text.x = element_text(face="bold")) +
scale_x_discrete(limits = c("Sex", "Species", "Strain"), expand = c(0.15, 0.05), position = "top") +
ggtitle("Species, strains and sex used (set2)")
freq_2 %>% filter(Freq != 0) %>% arrange(desc(Freq)) #collapesd table of values, without 0sAge at stress and enrichment
freq_ages1 <- as.data.frame(table(dat$Age_stress_exposure, dat$Age_EE_exposure)) %>% rename(Age_stress = Var1, Age_EE = Var2)
is_alluvia_form(as.data.frame(freq_ages1), axes = 1:2, silent = TRUE)
ggplot(data = freq_ages1,
aes(axis1 = Age_stress, axis2 = Age_EE, y = Freq)) +
geom_alluvium(aes(fill = Age_EE)) +
geom_stratum(aes(fill = Age_stress)) +
geom_text(stat = "stratum", aes(label = after_stat(stratum))) +
#theme_minimal() +
theme_void() +
theme(legend.position = "none",
plot.title = element_text(hjust = 0.5, vjust = 3),
axis.title.x = element_text(),
axis.text.x = element_text(face="bold")) +
scale_x_discrete(limits = c("Stress", "Enrichment"), expand = c(0.15, 0.05), position = "top") +
ggtitle("Life stages used (set1)")
Age at stress, enrichment, and assay
freq_ages2 <- as.data.frame(table(dat$Age_stress_exposure, dat$Age_EE_exposure, dat$Age_assay)) %>% rename(Age_stress = Var1, Age_EE = Var2, Age_assay = Var3)
is_alluvia_form(as.data.frame(freq_ages2), axes = 1:3, silent = TRUE)
ggplot(data = freq_ages2,
aes(axis1 = Age_stress, axis2 = Age_EE, axis3 = Age_assay, y = Freq)) +
geom_alluvium(aes(fill = Age_stress)) +
geom_flow() +
geom_stratum(aes(fill = Age_stress)) +
geom_text(stat = "stratum", aes(label = after_stat(stratum))) +
#theme_minimal() +
theme_void() +
theme(legend.position = "none",
plot.title = element_text(hjust = 0.5, vjust = 3),
axis.title.x = element_text(),
axis.text.x = element_text(face="bold")) +
scale_x_discrete(limits = c("Stress", "Enrichment", "Assay"), expand = c(0.15, 0.05), position = "top") +
ggtitle("Life stages used (set2)")
freq_ages2 %>% filter(Freq != 0) %>% arrange(desc(Freq)) #collapesd table of values, without 0sAge at stress, enrichment, assay, and exposure order
freq_ages3 <- as.data.frame(table(dat$Age_stress_exposure, dat$Age_EE_exposure, dat$Age_assay, dat$Exposure_order)) %>% rename(Age_stress = Var1, Age_EE = Var2, Age_assay = Var3, Order = Var4)
is_alluvia_form(as.data.frame(freq_ages3), axes = 1:4, silent = TRUE)
ggplot(data = freq_ages3,
aes(axis1 = Age_stress, axis2 = Age_EE, axis3 = Age_assay, axis4 = Order, y = Freq)) +
geom_alluvium(aes(fill = Age_stress)) +
geom_flow() +
geom_stratum(aes(fill = Age_stress)) +
geom_text(stat = "stratum", aes(label = after_stat(stratum))) +
#theme_minimal() +
theme_void() +
theme(legend.position = "none",
plot.title = element_text(hjust = 0.5, vjust = 3),
axis.title.x = element_text(),
axis.text.x = element_text(face="bold")) +
scale_x_discrete(limits = c("Stress", "Enrichment", "Assay", "Order"), expand = c(0.15, 0.05), position = "top") +
ggtitle("Life stages used and order (set3)")
freq_ages2 %>% filter(Freq != 0) %>% arrange(desc(Freq)) #collapesd table of values, without 0s
freq_ages3 %>% filter(Freq != 0) %>% arrange(desc(Freq)) #collapesd table of values, without 0sStress duration and type of stress
freq_stress1 <- as.data.frame(table(dat$Stress_duration, dat$Type_stress_exposure)) %>% rename(Stress_duration = Var1, Stress_type = Var2)
is_alluvia_form(as.data.frame(freq_stress1), axes = 1:2, silent = TRUE)
ggplot(data = freq_stress1,
aes(axis1 = Stress_duration, axis2 = Stress_type, y = Freq)) +
geom_alluvium(aes(fill = Stress_duration)) +
geom_stratum(aes(fill = Stress_duration)) +
geom_text(stat = "stratum", aes(label = after_stat(stratum))) +
#theme_minimal() +
theme_void() +
theme(legend.position = "none",
plot.title = element_text(hjust = 0.5, vjust = 3),
axis.title.x = element_text(),
axis.text.x = element_text(face="bold")) +
scale_x_discrete(limits = c("Duration", "Type"), expand = c(0.15, 0.05), position = "top") +
ggtitle("Stress exposure variables (set1)")
Age at stress, type of stress, stress duration
freq_stress2 <- as.data.frame(table(dat$Age_stress_exposure, dat$Stress_duration, dat$Type_stress_exposure)) %>% rename(Age_stress = Var1, Duration_stress = Var2, Type_stress = Var3)
is_alluvia_form(as.data.frame(freq_stress2), axes = 1:3, silent = TRUE)
ggplot(data = freq_stress2,
aes(axis1 = Age_stress, axis2 = Duration_stress, axis3 = Type_stress, y = Freq)) +
geom_alluvium(aes(fill = Age_stress)) +
geom_flow() +
geom_stratum(aes(fill = Age_stress)) +
geom_text(stat = "stratum", aes(label = after_stat(stratum))) +
#theme_minimal() +
theme_void() +
theme(legend.position = "none",
plot.title = element_text(hjust = 0.5, vjust = 3),
axis.title.x = element_text(),
axis.text.x = element_text(face="bold")) +
scale_x_discrete(limits = c("Age", "Duration", "Type"), expand = c(0.1, 0.1), position = "top") +
ggtitle("Stress exposure variables (set2)")
Exercise and social enrichment
freq_EE1 <- as.data.frame(table(dat$EE_exercise, dat$EE_social)) %>% rename(EE_exercise = Var1, EE_social = Var2)
is_alluvia_form(as.data.frame(freq_stress1), axes = 1:2, silent = TRUE)
ggplot(data = freq_EE1,
aes(axis1 = EE_exercise, axis2 = EE_social, y = Freq)) +
geom_alluvium(aes(fill = EE_exercise)) +
geom_stratum(aes(fill = EE_exercise)) +
geom_text(stat = "stratum", aes(label = after_stat(stratum))) +
#theme_minimal() +
theme_void() +
theme(legend.position = "none",
plot.title = element_text(hjust = 0.5, vjust = 3),
axis.title.x = element_text(),
axis.text.x = element_text(face="bold")) +
scale_x_discrete(limits = c("Exercise", "Social"), expand = c(0.15, 0.05), position = "top") +
ggtitle("EE exposure variables (set1)")
Age at enrichment, exercise and social enrichment
freq_EE2 <- as.data.frame(table(dat$Age_EE_exposure, dat$EE_exercise, dat$EE_social)) %>% rename(Age_EE = Var1, EE_exercise = Var2, EE_social = Var3)
is_alluvia_form(as.data.frame(freq_EE2), axes = 1:3, silent = TRUE)
ggplot(data = freq_EE2,
aes(axis1 = Age_EE, axis2 = EE_exercise, axis3 = EE_social, y = Freq)) +
geom_alluvium(aes(fill = Age_EE)) +
geom_flow() +
geom_stratum(aes(fill = Age_EE)) +
geom_text(stat = "stratum", aes(label = after_stat(stratum))) +
#theme_minimal() +
theme_void() +
theme(legend.position = "none",
plot.title = element_text(hjust = 0.5, vjust = 3),
axis.title.x = element_text(),
axis.text.x = element_text(face="bold")) +
scale_x_discrete(limits = c("Age", "Exercise", "Social"), expand = c(0.1, 0.1), position = "top") +
ggtitle("EE exposure variables (set2)")
Learning vs memory and type of reinforcement
freq_assay1 <- as.data.frame(table(dat$Learning_vs_memory, dat$Type_reinforcement)) %>% rename(Learning_Memory = Var1, Reinforcement = Var2)
is_alluvia_form(as.data.frame(freq_assay1), axes = 1:2, silent = TRUE)
ggplot(data = freq_assay1,
aes(axis1 = Learning_Memory, axis2 = Reinforcement, y = Freq)) +
geom_alluvium(aes(fill = Learning_Memory)) +
geom_stratum(aes(fill = Learning_Memory)) +
geom_text(stat = "stratum", aes(label = after_stat(stratum))) +
#theme_minimal() +
theme_void() +
theme(legend.position = "none",
plot.title = element_text(hjust = 0.5, vjust = 3),
axis.title.x = element_text(),
axis.text.x = element_text(face="bold")) +
scale_x_discrete(limits = c("Learning_Memory", "Reinforcement"), expand = c(0.15, 0.05), position = "top") +
ggtitle("Behavioural assay variables (set1)")
Age at assay, learning vs memory, and type of reinforcement
freq_assay2 <- as.data.frame(table(dat$Age_assay, dat$Learning_vs_memory, dat$Type_reinforcement)) %>% rename(Age_assay = Var1, Learning_Memory = Var2, Reinforcement = Var3)
is_alluvia_form(as.data.frame(freq_assay2), axes = 1:3, silent = TRUE)
ggplot(data = freq_assay2,
aes(axis1 = Age_assay, axis2 = Learning_Memory, axis3 = Reinforcement, y = Freq)) +
geom_alluvium(aes(fill = Age_assay)) +
geom_flow() +
geom_stratum(aes(fill = Age_assay)) +
geom_text(stat = "stratum", aes(label = after_stat(stratum))) +
#theme_minimal() +
theme_void() +
theme(legend.position = "none",
plot.title = element_text(hjust = 0.5, vjust = 3),
axis.title.x = element_text(),
axis.text.x = element_text(face="bold")) +
scale_x_discrete(limits = c("Age", "Learning_Memory", "Reinforcement"), expand = c(0.1, 0.1), position = "top") +
ggtitle("Behavioural assay variables (set2)")
Type of learning, learning vs memory, type of reinforcement
freq_assay3 <- as.data.frame(table(dat$Type_learning, dat$Learning_vs_memory, dat$Type_reinforcement)) %>% rename(Type = Var1, Learning_Memory = Var2, Reinforcement = Var3)
is_alluvia_form(as.data.frame(freq_assay3), axes = 1:3, silent = TRUE)
ggplot(data = freq_assay3,
aes(axis1 = Type, axis2 = Learning_Memory, axis3 = Reinforcement, y = Freq)) +
geom_alluvium(aes(fill = Type)) +
geom_flow() +
geom_stratum(aes(fill = Type)) +
geom_text(stat = "stratum", aes(label = after_stat(stratum))) +
#theme_minimal() +
theme_void() +
theme(legend.position = "none",
plot.title = element_text(hjust = 0.5, vjust = 3),
axis.title.x = element_text(),
axis.text.x = element_text(face="bold")) +
scale_x_discrete(limits = c("Type", "Learning_Memory", "Reinforcement"), expand = c(0.1, 0.1), position = "top") +
ggtitle("Behavioural assay variables (set3)")
Modelling with lnRR
Environmental enrichment
Intercept model
Learning and memory are significantly improved when housed with environmnetal enrichment
mod_E0 <- rma.mv(yi = lnRR_Ea, V = VCV_E, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat)
summary(mod_E0) ##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -20.3125 40.6249 48.6249 58.6683 49.0900
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0044 0.0665 30 no Study_ID
## sigma^2.2 0.0485 0.2203 92 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Heterogeneity:
## Q(df = 91) = 897.6814, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## 0.2146 0.0339 6.3344 91 <.0001 0.1473 0.2818 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1i2_ml(mod_E0) ## I2_total I2_Study_ID I2_ES_ID I2_Strain
## 9.398214e-01 7.853525e-02 8.612861e-01 7.246533e-10orchard_plot(mod_E0, mod = "Int", xlab = "lnRR", alpha=0.4) + # Orchard plot
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5)+ # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2)+ # confidence intervals
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_colour_manual(values = "darkorange")+ # change colours
scale_fill_manual(values="darkorange")+
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 24), # change font sizes
legend.title = element_text(size = 15),
legend.text = element_text(size = 13)) 
Metaregression
Uni-Moderator metaregression
Type of learning
The type of learning/memory response
dat1 <- filter(dat, Type_learning %in% c("Recognition", "Habituation", "Conditioning"))
VCV_E1 <- impute_covariance_matrix(vi = dat1$lnRRV_E, cluster = dat$Study_ID, r = 0.5)
mod_E1 <- rma.mv(yi = lnRR_Ea, V = VCV_E1, mod = ~Type_learning-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat1)
summary(mod_E1)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -18.0617 36.1234 48.1234 63.0553 49.1478
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0081 0.0898 30 no Study_ID
## sigma^2.2 0.0434 0.2083 92 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 89) = 790.2293, p-val < .0001
##
## Test of Moderators (coefficients 1:3):
## F(df1 = 3, df2 = 89) = 14.6747, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb
## Type_learningConditioning 0.2514 0.0383 6.5700 89 <.0001 0.1754
## Type_learningHabituation 0.2267 0.1216 1.8639 89 0.0656 -0.0150
## Type_learningRecognition 0.0635 0.0706 0.8997 89 0.3707 -0.0768
## ci.ub
## Type_learningConditioning 0.3275 ***
## Type_learningHabituation 0.4684 .
## Type_learningRecognition 0.2039
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1r2_ml(mod_E1) ## R2_marginal R2_coditional
## 0.08104773 1.00000000Learning_E <- orchard_plot(mod_E1, mod = "Type_learning", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
xlim(-0.5, 2) +
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 15), # change font sizes
legend.title = element_text(size = 10),
legend.text = element_text(size = 10))
Learning_E
Learning vs Memory
Is the assay broadly measuring learning or memory?
mod_E2 <- rma.mv(yi = lnRR_Ea, V = VCV_E, mod = ~Learning_vs_memory-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat)
summary(mod_E2) ##
## Multivariate Meta-Analysis Model (k = 85; method: REML)
##
## logLik Deviance AIC BIC AICc
## -11.6265 23.2529 33.2529 45.3471 34.0322
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0077 0.0877 30 no Study_ID
## sigma^2.2 0.0310 0.1760 85 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 83) = 652.9299, p-val < .0001
##
## Test of Moderators (coefficients 1:2):
## F(df1 = 2, df2 = 83) = 21.0340, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb
## Learning_vs_memoryLearning 0.2458 0.0475 5.1789 83 <.0001 0.1514
## Learning_vs_memoryMemory 0.1884 0.0360 5.2368 83 <.0001 0.1169
## ci.ub
## Learning_vs_memoryLearning 0.3402 ***
## Learning_vs_memoryMemory 0.2600 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1r2_ml(mod_E2) ## R2_marginal R2_coditional
## 0.01866131 1.00000000LvsM_E<- orchard_plot(mod_E2, mod = "Learning_vs_memory", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
xlim(-0.5, 2) +
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 15), # change font sizes
legend.title = element_text(size = 10),
legend.text = element_text(size = 10))
LvsM_E
Type of reinforcement
The type of cue used
dat2 <- filter(dat, Type_reinforcement %in% c("Appetitive", "Aversive", "Not applicable"))
VCV_E2 <- impute_covariance_matrix(vi = dat2$lnRRV_E, cluster = dat2$Study_ID, r = 0.5)
mod_E3 <- rma.mv(yi = lnRR_Ea, V = VCV_E2, mod = ~ Type_reinforcement-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat2)
summary(mod_E3)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -18.1008 36.2017 48.2017 63.1335 49.2261
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0061 0.0778 30 no Study_ID
## sigma^2.2 0.0449 0.2119 92 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 89) = 780.4926, p-val < .0001
##
## Test of Moderators (coefficients 1:3):
## F(df1 = 3, df2 = 89) = 15.1398, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb
## Type_reinforcementAppetitive 0.1949 0.0721 2.7023 89 0.0082 0.0516
## Type_reinforcementAversive 0.2704 0.0439 6.1559 89 <.0001 0.1831
## Type_reinforcementNot applicable 0.1082 0.0600 1.8030 89 0.0748 -0.0110
## ci.ub
## Type_reinforcementAppetitive 0.3382 **
## Type_reinforcementAversive 0.3577 ***
## Type_reinforcementNot applicable 0.2274 .
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1r2_ml(mod_E3) ## R2_marginal R2_coditional
## 0.08421309 1.00000000Reinforcement_E <-orchard_plot(mod_E3, mod = "Type_reinforcement", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
xlim(-0.5, 2) +
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 15), # change font sizes
legend.title = element_text(size = 10),
legend.text = element_text(size = 10))
Reinforcement_E
Exercise enrichment
Does the form of enrichment include exercise through a running wheel or treadmill?
mod_E5<- rma.mv(yi = lnRR_Ea, V = VCV_E, mod = ~EE_exercise-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat)
summary(mod_E5)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -20.4952 40.9905 50.9905 63.4895 51.7048
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0054 0.0734 30 no Study_ID
## sigma^2.2 0.0487 0.2208 92 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 90) = 867.0038, p-val < .0001
##
## Test of Moderators (coefficients 1:2):
## F(df1 = 2, df2 = 90) = 19.4577, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## EE_exerciseExercise 0.2156 0.0421 5.1141 90 <.0001 0.1318 0.2993
## EE_exerciseNo exercise 0.2146 0.0601 3.5723 90 0.0006 0.0952 0.3339
##
## EE_exerciseExercise ***
## EE_exerciseNo exercise ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1r2_ml(mod_E5) ## R2_marginal R2_coditional
## 4.108541e-06 1.000000e+00Exercise_E <-orchard_plot(mod_E5, mod = "EE_exercise", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
xlim(-0.5, 2) +
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 15), # change font sizes
legend.title = element_text(size = 10),
legend.text = element_text(size = 10))
Exercise_E
Age of enrichment
The age at which the individuals were exposed to environmental enrichment.
dat6 <- filter(dat, Age_EE_exposure %in% c("Adult", "Adolescent"))
VCV_E6 <- impute_covariance_matrix(vi = dat6$lnRRV_E, cluster = dat6$Study_ID, r = 0.5)
mod_E6 <- rma.mv(yi = lnRR_Ea, V = VCV_E6, mod = ~Age_EE_exposure-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat6)
summary(mod_E6) ##
## Multivariate Meta-Analysis Model (k = 88; method: REML)
##
## logLik Deviance AIC BIC AICc
## -17.2225 34.4450 44.4450 56.7168 45.1950
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0047 0.0682 29 no Study_ID
## sigma^2.2 0.0457 0.2137 88 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 86) = 834.1449, p-val < .0001
##
## Test of Moderators (coefficients 1:2):
## F(df1 = 2, df2 = 86) = 21.8297, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## Age_EE_exposureAdolescent 0.2115 0.0389 5.4335 86 <.0001 0.1341 0.2889
## Age_EE_exposureAdult 0.2572 0.0684 3.7599 86 0.0003 0.1212 0.3932
##
## Age_EE_exposureAdolescent ***
## Age_EE_exposureAdult ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1r2_ml(mod_E6) ## R2_marginal R2_coditional
## 0.006503496 0.999999999Age_E<- orchard_plot(mod_E6, mod = "Age_EE_exposure", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
xlim(-0.5, 2) +
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 15), # change font sizes
legend.title = element_text(size = 10),
legend.text = element_text(size = 10))
Age_E
multi-moderator Full model
# filter data so that all K < 5 are removed
dat_Efm <- dat %>%
filter(Type_learning %in% c("Recognition", "Habituation", "Conditioning"),
Type_reinforcement %in% c("Appetitive", "Aversive", "Not applicable"),
EE_social %in% c("Social", "Non-social"),
Age_EE_exposure %in% c("Adult", "Adolescent"))
VCV_Efm <- impute_covariance_matrix(vi = dat_Efm$lnRRV_E, cluster = dat$Study_ID, r = 0.5)
mod_Efm <- rma.mv(yi = lnRR_Sa, V = VCV_Efm, mod = ~Type_learning-1 + Learning_vs_memory-1 + Type_reinforcement-1 + EE_social-1 + EE_exercise-1 + Age_EE_exposure-1 , random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat_Efm)
#summary(mod_Efm)
#r2_ml(mod_Efm)
res_Efm <- dredge(mod_Efm, trace=2)##
|
| | 0%
|
|= | 2%
|
|== | 3%
|
|=== | 5%
|
|==== | 6%
|
|===== | 8%
|
|======= | 9%
|
|======== | 11%
|
|========= | 12%
|
|========== | 14%
|
|=========== | 16%
|
|============ | 17%
|
|============= | 19%
|
|============== | 20%
|
|=============== | 22%
|
|================ | 23%
|
|================== | 25%
|
|=================== | 27%
|
|==================== | 28%
|
|===================== | 30%
|
|====================== | 31%
|
|======================= | 33%
|
|======================== | 34%
|
|========================= | 36%
|
|========================== | 38%
|
|=========================== | 39%
|
|============================ | 41%
|
|============================== | 42%
|
|=============================== | 44%
|
|================================ | 45%
|
|================================= | 47%
|
|================================== | 48%
|
|=================================== | 50%
|
|==================================== | 52%
|
|===================================== | 53%
|
|====================================== | 55%
|
|======================================= | 56%
|
|======================================== | 58%
|
|========================================== | 59%
|
|=========================================== | 61%
|
|============================================ | 62%
|
|============================================= | 64%
|
|============================================== | 66%
|
|=============================================== | 67%
|
|================================================ | 69%
|
|================================================= | 70%
|
|================================================== | 72%
|
|=================================================== | 73%
|
|==================================================== | 75%
|
|====================================================== | 77%
|
|======================================================= | 78%
|
|======================================================== | 80%
|
|========================================================= | 81%
|
|========================================================== | 83%
|
|=========================================================== | 84%
|
|============================================================ | 86%
|
|============================================================= | 88%
|
|============================================================== | 89%
|
|=============================================================== | 91%
|
|================================================================= | 92%
|
|================================================================== | 94%
|
|=================================================================== | 95%
|
|==================================================================== | 97%
|
|===================================================================== | 98%res_Efm2<- subset(res_Efm, delta <= 6, recalc.weights=FALSE)
importance(res_Efm2)## Learning_vs_memory Age_EE_exposure EE_exercise
## Sum of weights: 0.94 0.52 0.24
## N containing models: 20 10 8
## Type_reinforcement EE_social Type_learning
## Sum of weights: 0.22 0.21 0.14
## N containing models: 7 7 6Publication bias & sensitivity analysis
Funnel_E<-funnel(mod_Efm, xlab = "lnRR", ylab = "Standard Error")
#year published was scaled previously under stress PB
dat_Efm$sqrt_inv_e_n <- with(dat_Efm, sqrt(1/CC_n + 1/EC_n + 1/ES_n + 1/CS_n))
PB_MR_E<- rma.mv(lnRR_Sa, lnRRV_S, mods = ~1 + sqrt_inv_e_n + Year_published + Type_learning + Type_reinforcement + EE_social + EE_exercise + Age_stress_exposure, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain), method = "REML", test = "t",
data = dat_Efm)
estimates_PB_MR_E<- estimates.CI(PB_MR_E)
estimates_PB_MR_E| estimate | mean | lower | upper |
|---|---|---|---|
| intrcpt | -13.9147389 | -59.3987797 | 31.5693018 |
| sqrt_inv_e_n | -0.1711961 | -1.0441687 | 0.7017764 |
| Year_published | 0.0068520 | -0.0157749 | 0.0294790 |
| Type_learningHabituation | -0.5828922 | -1.0088024 | -0.1569819 |
| Type_learningRecognition | -0.1068229 | -0.4651361 | 0.2514902 |
| Type_reinforcementAversive | 0.2077553 | -0.1500255 | 0.5655360 |
| Type_reinforcementNot applicable | 0.3274613 | -0.1646052 | 0.8195278 |
| EE_socialSocial | 0.0540916 | -0.1861944 | 0.2943777 |
| EE_exerciseNo exercise | -0.0193952 | -0.2891022 | 0.2503118 |
| Age_stress_exposureAdult | -0.1438520 | -0.5799393 | 0.2922354 |
| Age_stress_exposureEarly postnatal | -0.0499447 | -0.4906291 | 0.3907397 |
| Age_stress_exposurePrenatal | -0.1335549 | -0.5816002 | 0.3144903 |
#no effect of inv_effective sample size or year published
#leave-one-out analysis
dat$Study_ID <- as.factor(dat$Study_ID)
LeaveOneOut_effectsize <- list()
for(i in 1:length(levels(dat$Study_ID))){
LeaveOneOut_effectsize[[i]] <- rma.mv(yi = lnRR_Ea, V = lnRRV_E,
random = list(~1 | Study_ID,~1| ES_ID, ~1 | Strain),
method = "REML", data = dat[dat$Study_ID
!= levels(dat$Study_ID)[i], ])}
# writing function for extracting est, ci.lb, and ci.ub from all models
est.func <- function(mod_E0){
df <- data.frame(est = mod_E0$b, lower = mod_E0$ci.lb, upper = mod_E0$ci.ub)
return(df)
}
#using dplyr to form data frame
MA_CVR <- lapply(LeaveOneOut_effectsize, function(x) est.func(x))%>% bind_rows %>% mutate(left_out = levels(dat$Study_ID))
#telling ggplot to stop reordering factors
MA_CVR$left_out<- as.factor(MA_CVR$left_out)
MA_CVR$left_out<-factor(MA_CVR$left_out, levels = MA_CVR$left_out)
#plotting
leaveoneout_E <- ggplot(MA_CVR) +
geom_hline(yintercept = 0, lty = 2, lwd = 1) +
geom_hline(yintercept = mod_E0$ci.lb, lty = 3, lwd = 0.75, colour = "black") +
geom_hline(yintercept = mod_E0$b, lty = 1, lwd = 0.75, colour = "black") +
geom_hline(yintercept = mod_E0$ci.ub, lty = 3, lwd = 0.75, colour = "black") +
geom_pointrange(aes(x = left_out, y = est, ymin = lower, ymax = upper)) +
xlab("Study left out") +
ylab("lnRR, 95% CI") +
coord_flip() +
theme(panel.grid.minor = element_blank())+
theme_bw() + theme(panel.grid.major = element_blank()) +
theme(panel.grid.minor.x = element_blank() ) +
theme(axis.text.y = element_text(size = 6))
leaveoneout_E
dat$Study_ID <- as.integer(dat$Study_ID)Stress
Intercept model
Learning and memory are significantly reduced due to stress. High heterogeneity
mod_S0 <- rma.mv(yi = lnRR_Sa, V = VCV_S, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat)
summary(mod_S0) ##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -22.6902 45.3804 53.3804 63.4239 53.8456
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0107 0.1036 30 no Study_ID
## sigma^2.2 0.0504 0.2245 92 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Heterogeneity:
## Q(df = 91) = 933.9737, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## -0.1155 0.0378 -3.0591 91 0.0029 -0.1906 -0.0405 **
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1i2_ml(mod_S0) ## I2_total I2_Study_ID I2_ES_ID I2_Strain
## 9.474395e-01 1.662608e-01 7.811787e-01 2.947034e-10orchard_plot(mod_S0, mod = "Int", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 24), # change font sizes
legend.title = element_text(size = 15),
legend.text = element_text(size = 13)) 
###Metaregression
Uni-moderator metaregression
Type of learning
The type of learning/memory response
dat$Type_learning<-as.factor(dat$Type_learning)
VCV_S1 <- impute_covariance_matrix(vi = dat1$lnRRV_S, cluster = dat$Study_ID, r = 0.5)
mod_S1 <- rma.mv(yi = lnRR_Sa, V = VCV_S1, mod = ~Type_learning-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat1)
summary(mod_S1)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -17.3159 34.6317 46.6317 61.5635 47.6561
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0203 0.1423 30 no Study_ID
## sigma^2.2 0.0351 0.1875 92 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 89) = 732.1642, p-val < .0001
##
## Test of Moderators (coefficients 1:3):
## F(df1 = 3, df2 = 89) = 7.3668, p-val = 0.0002
##
## Model Results:
##
## estimate se tval df pval ci.lb
## Type_learningConditioning -0.1052 0.0425 -2.4791 89 0.0151 -0.1896
## Type_learningHabituation -0.5273 0.1191 -4.4285 89 <.0001 -0.7639
## Type_learningRecognition -0.0490 0.0718 -0.6819 89 0.4971 -0.1917
## ci.ub
## Type_learningConditioning -0.0209 *
## Type_learningHabituation -0.2907 ***
## Type_learningRecognition 0.0937
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1r2_ml(mod_S1) ## R2_marginal R2_coditional
## 0.1974766 1.0000000Learning_S <-orchard_plot(mod_S1, mod = "Type_learning", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 15), # change font sizes
legend.title = element_text(size = 10),
legend.text = element_text(size = 10))
Learning_S
Learning vs Memory
Is the assay broadly measuring learning or memory?
mod_S2 <- rma.mv(yi = lnRR_Sa, V = VCV_S, mod = ~Learning_vs_memory-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat)
summary(mod_S2) ##
## Multivariate Meta-Analysis Model (k = 85; method: REML)
##
## logLik Deviance AIC BIC AICc
## -11.8032 23.6065 33.6065 45.7007 34.3857
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0261 0.1615 30 no Study_ID
## sigma^2.2 0.0267 0.1635 85 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 83) = 676.7747, p-val < .0001
##
## Test of Moderators (coefficients 1:2):
## F(df1 = 2, df2 = 83) = 2.9221, p-val = 0.0594
##
## Model Results:
##
## estimate se tval df pval ci.lb
## Learning_vs_memoryLearning -0.0573 0.0532 -1.0785 83 0.2840 -0.1631
## Learning_vs_memoryMemory -0.1053 0.0436 -2.4135 83 0.0180 -0.1920
## ci.ub
## Learning_vs_memoryLearning 0.0484
## Learning_vs_memoryMemory -0.0185 *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1r2_ml(mod_S2) ## R2_marginal R2_coditional
## 0.009625517 1.000000000LvsM_S <- orchard_plot(mod_S2, mod = "Learning_vs_memory", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 15), # change font sizes
legend.title = element_text(size = 10),
legend.text = element_text(size = 10))
LvsM_S 
Type of reinforcement
The type of cue used
VCV_S2 <- VCV_E <- impute_covariance_matrix(vi = dat2$lnRRV_S, cluster = dat$Study_ID, r = 0.5)
mod_S3 <- rma.mv(yi = lnRR_Sa, V = VCV_S2, mod = ~ Type_reinforcement-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat2)
summary(mod_S3)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -22.1655 44.3311 56.3311 71.2629 57.3555
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0126 0.1124 30 no Study_ID
## sigma^2.2 0.0506 0.2248 92 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 89) = 911.1858, p-val < .0001
##
## Test of Moderators (coefficients 1:3):
## F(df1 = 3, df2 = 89) = 3.6161, p-val = 0.0162
##
## Model Results:
##
## estimate se tval df pval
## Type_reinforcementAppetitive -0.1980 0.0836 -2.3685 89 0.0200
## Type_reinforcementAversive -0.0747 0.0489 -1.5264 89 0.1305
## Type_reinforcementNot applicable -0.1386 0.0660 -2.0995 89 0.0386
## ci.lb ci.ub
## Type_reinforcementAppetitive -0.3640 -0.0319 *
## Type_reinforcementAversive -0.1719 0.0225
## Type_reinforcementNot applicable -0.2698 -0.0074 *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1r2_ml(mod_S3) ## R2_marginal R2_coditional
## 0.03676953 1.00000000Reinforcement_S <-orchard_plot(mod_S3, mod = "Type_reinforcement", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 15), # change font sizes
legend.title = element_text(size = 10),
legend.text = element_text(size = 10))
Reinforcement_S
Type of stress
The type of stress manipulation
dat3 <- filter(dat, Type_stress_exposure %in% c("Restraint", "Noise", "MS", "Combination"))
VCV_S3 <- impute_covariance_matrix(vi = dat3$lnRRV_S, cluster = dat$Study_ID, r = 0.5)
mod_S4 <- rma.mv(yi = lnRR_Sa, V = VCV_S3, mod = ~Type_stress_exposure-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat3)
summary(mod_S4) ##
## Multivariate Meta-Analysis Model (k = 85; method: REML)
##
## logLik Deviance AIC BIC AICc
## -21.3621 42.7241 56.7241 73.4853 58.2584
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0236 0.1537 25 no Study_ID
## sigma^2.2 0.0527 0.2295 85 no ES_ID
## sigma^2.3 0.0000 0.0000 4 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 81) = 1030.5305, p-val < .0001
##
## Test of Moderators (coefficients 1:4):
## F(df1 = 4, df2 = 81) = 2.0257, p-val = 0.0985
##
## Model Results:
##
## estimate se tval df pval ci.lb
## Type_stress_exposureCombination -0.0558 0.1132 -0.4927 81 0.6236 -0.2811
## Type_stress_exposureMS -0.0607 0.0691 -0.8780 81 0.3825 -0.1982
## Type_stress_exposureNoise -0.1468 0.1248 -1.1765 81 0.2428 -0.3950
## Type_stress_exposureRestraint -0.2175 0.0878 -2.4764 81 0.0154 -0.3922
## ci.ub
## Type_stress_exposureCombination 0.1695
## Type_stress_exposureMS 0.0768
## Type_stress_exposureNoise 0.1015
## Type_stress_exposureRestraint -0.0427 *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1r2_ml(mod_S4)## R2_marginal R2_coditional
## 0.05175205 1.00000000Stressor<- orchard_plot(mod_S4, mod = "Type_stress_exposure", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 15), # change font sizes
legend.title = element_text(size = 10),
legend.text = element_text(size = 10))
Stressor
Age of stress
VCV_S3a <- impute_covariance_matrix(vi = dat$lnRRV_S, cluster = dat$Study_ID, r = 0.5)
mod_S5 <-rma.mv(yi = lnRR_Sa, V = VCV_S3a, mod = ~Age_stress_exposure-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat)
summary(mod_S5) ##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -20.8850 41.7699 55.7699 73.1113 57.1699
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0021 0.0455 30 no Study_ID
## sigma^2.2 0.0531 0.2304 92 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 88) = 830.4810, p-val < .0001
##
## Test of Moderators (coefficients 1:4):
## F(df1 = 4, df2 = 88) = 4.4779, p-val = 0.0024
##
## Model Results:
##
## estimate se tval df pval
## Age_stress_exposureAdolescent 0.0244 0.1336 0.1825 88 0.8556
## Age_stress_exposureAdult -0.2481 0.0693 -3.5790 88 0.0006
## Age_stress_exposureEarly postnatal -0.0645 0.0461 -1.3997 88 0.1651
## Age_stress_exposurePrenatal -0.1379 0.0782 -1.7634 88 0.0813
## ci.lb ci.ub
## Age_stress_exposureAdolescent -0.2411 0.2899
## Age_stress_exposureAdult -0.3859 -0.1104 ***
## Age_stress_exposureEarly postnatal -0.1561 0.0271
## Age_stress_exposurePrenatal -0.2932 0.0175 .
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1r2_ml(mod_S5) ## R2_marginal R2_coditional
## 0.0971388 1.0000000Age_S <- orchard_plot(mod_S5, mod = "Age_stress_exposure", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 15), # change font sizes
legend.title = element_text(size = 10),
legend.text = element_text(size = 10))
Age_S 
Stess duration
How long was the stress applied for (chronic = every day for 7 days or more)? This has the highest marginal R2
dat4 <- filter(dat, Stress_duration %in% c("Chronic", "Acute"))
VCV_S4 <- impute_covariance_matrix(vi = dat4$lnRRV_S, cluster = dat$Study_ID, r = 0.5)
mod_S6 <-rma.mv(yi = lnRR_Sa, V = VCV_S4, mod = ~Stress_duration-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat4)
summary(mod_S6) ##
## Multivariate Meta-Analysis Model (k = 89; method: REML)
##
## logLik Deviance AIC BIC AICc
## -23.6341 47.2682 57.2682 69.5978 58.0090
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0173 0.1314 29 no Study_ID
## sigma^2.2 0.0514 0.2267 89 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 87) = 1164.1990, p-val < .0001
##
## Test of Moderators (coefficients 1:2):
## F(df1 = 2, df2 = 87) = 4.9979, p-val = 0.0088
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## Stress_durationAcute 0.0025 0.0866 0.0292 87 0.9768 -0.1695 0.1746
## Stress_durationChronic -0.1509 0.0477 -3.1593 87 0.0022 -0.2458 -0.0559
##
## Stress_durationAcute
## Stress_durationChronic **
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1r2_ml(mod_S6) ## R2_marginal R2_coditional
## 0.05879435 1.00000000Duration_S <- orchard_plot(mod_S6, mod = "Stress_duration", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 15), # change font sizes
legend.title = element_text(size = 10),
legend.text = element_text(size = 10))
Duration_S 
Multi-moderator fullmodel
#selecting moderator levels with k >=5
dat_Sfm <- dat %>%
filter(Type_learning %in% c("Recognition", "Habituation", "Conditioning"),
Type_reinforcement %in% c("Appetitive", "Aversive", "Not applicable"),
Type_stress_exposure %in% c("Restraint", "Noise", "MS", "Combination"),
Stress_duration %in% c("Chronic", "Acute"))
VCV_Sfm <- impute_covariance_matrix(vi = dat_Sfm$lnRRV_E, cluster = dat$Study_ID, r = 0.5)
mod_Sfm <- rma.mv(yi = lnRR_Sa, V = VCV_Sfm, mod = ~ Type_learning -1 + Learning_vs_memory + Type_reinforcement + Type_stress_exposure + Age_stress_exposure + Stress_duration, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat_Sfm)
#summary(mod_Sfm)
#r2_ml(mod_Sfm)
res_Sfm <- dredge(mod_Sfm, trace=2)##
|
| | 0%
|
|= | 2%
|
|== | 3%
|
|=== | 5%
|
|==== | 6%
|
|===== | 8%
|
|======= | 9%
|
|======== | 11%
|
|========= | 12%
|
|========== | 14%
|
|=========== | 16%
|
|============ | 17%
|
|============= | 19%
|
|============== | 20%
|
|=============== | 22%
|
|================ | 23%
|
|================== | 25%
|
|=================== | 27%
|
|==================== | 28%
|
|===================== | 30%
|
|====================== | 31%
|
|======================= | 33%
|
|======================== | 34%
|
|========================= | 36%
|
|========================== | 38%
|
|=========================== | 39%
|
|============================ | 41%
|
|============================== | 42%
|
|=============================== | 44%
|
|================================ | 45%
|
|================================= | 47%
|
|================================== | 48%
|
|=================================== | 50%
|
|==================================== | 52%
|
|===================================== | 53%
|
|====================================== | 55%
|
|======================================= | 56%
|
|======================================== | 58%
|
|========================================== | 59%
|
|=========================================== | 61%
|
|============================================ | 62%
|
|============================================= | 64%
|
|============================================== | 66%
|
|=============================================== | 67%
|
|================================================ | 69%
|
|================================================= | 70%
|
|================================================== | 72%
|
|=================================================== | 73%
|
|==================================================== | 75%
|
|====================================================== | 77%
|
|======================================================= | 78%
|
|======================================================== | 80%
|
|========================================================= | 81%
|
|========================================================== | 83%
|
|=========================================================== | 84%
|
|============================================================ | 86%
|
|============================================================= | 88%
|
|============================================================== | 89%
|
|=============================================================== | 91%
|
|================================================================= | 92%
|
|================================================================== | 94%
|
|=================================================================== | 95%
|
|==================================================================== | 97%
|
|===================================================================== | 98%res_Sfm2<- subset(res_Sfm, delta <= 6, recalc.weights=FALSE)
importance(res_Sfm2) ## Learning_vs_memory Stress_duration Type_learning
## Sum of weights: 0.95 0.95 0.20
## N containing models: 8 8 4
## Age_stress_exposure Type_stress_exposure
## Sum of weights: 0.16 0.14
## N containing models: 2 2
## Type_reinforcement
## Sum of weights: 0.14
## N containing models: 2Publication bias & sensitivity analysis
funnel(mod_Sfm, xlab = "lnRR", ylab = "Standard Error")
#calculating inv effective sample size for use in full meta-regression
dat_Sfm$sqrt_inv_e_n <- with(dat_Sfm, sqrt(1/CC_n + 1/EC_n + 1/ES_n + 1/CS_n))
#time lag bias and eggers regression
dat_Sfm$c_Year_published <- as.vector(scale(dat_Sfm$Year_published, scale = F))
PB_MR_S<- rma.mv(lnRR_Sa, lnRRV_S, mods = ~1 + sqrt_inv_e_n + c_Year_published + Type_learning + Type_reinforcement + Type_stress_exposure + Age_stress_exposure, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
method = "REML",
test = "t",
data = dat_Sfm,
control=list(optimizer="optim", optmethod="Nelder-Mead"))
estimates_PB_MR_S<- estimates.CI(PB_MR_S)
estimates_PB_MR_S| estimate | mean | lower | upper |
|---|---|---|---|
| intrcpt | -0.0534910 | -0.9384247 | 0.8314427 |
| sqrt_inv_e_n | 0.1476637 | -1.2155241 | 1.5108514 |
| c_Year_published | 0.0055548 | -0.0243577 | 0.0354673 |
| Type_learningHabituation | -0.6361107 | -1.0644764 | -0.2077451 |
| Type_learningRecognition | -0.1068719 | -0.4754180 | 0.2616742 |
| Type_reinforcementAversive | 0.1220611 | -0.1344910 | 0.3786133 |
| Type_reinforcementNot applicable | 0.2373222 | -0.1952699 | 0.6699142 |
| Type_stress_exposureMS | -0.0009972 | -0.8170353 | 0.8150410 |
| Type_stress_exposureNoise | 0.1814221 | -0.8280961 | 1.1909402 |
| Type_stress_exposureRestraint | -0.1312725 | -0.5964008 | 0.3338558 |
| Age_stress_exposureAdult | -0.2037280 | -0.6376774 | 0.2302214 |
| Age_stress_exposureEarly postnatal | -0.2129606 | -1.0867878 | 0.6608667 |
| Age_stress_exposurePrenatal | -0.1999017 | -0.7352586 | 0.3354553 |
#no effect of inv_effective sample size or year published
#leave-one-out analysis
dat$Study_ID <- as.factor(dat$Study_ID)
LeaveOneOut_effectsize <- list()
for(i in 1:length(levels(dat$Study_ID))){
LeaveOneOut_effectsize[[i]] <- rma.mv(yi = lnRR_Sa, V = lnRRV_S,
random = list(~1 | Study_ID,~1| ES_ID, ~1 | Strain),
method = "REML", data = dat[dat$Study_ID
!= levels(dat$Study_ID)[i], ])}
# writing function for extracting est, ci.lb, and ci.ub from all models
est.func <- function(mod_E0){
df <- data.frame(est = mod_E0$b, lower = mod_E0$ci.lb, upper = mod_E0$ci.ub)
return(df)
}
#using dplyr to form data frame
MA_CVR <- lapply(LeaveOneOut_effectsize, function(x) est.func(x))%>% bind_rows %>% mutate(left_out = levels(dat$Study_ID))
#telling ggplot to stop reordering factors
MA_CVR$left_out<- as.factor(MA_CVR$left_out)
MA_CVR$left_out<-factor(MA_CVR$left_out, levels = MA_CVR$left_out)
#plotting
leaveoneout_S <- ggplot(MA_CVR) +
geom_hline(yintercept = 0, lty = 2, lwd = 1) +
geom_hline(yintercept = mod_S0$ci.lb, lty = 3, lwd = 0.75, colour = "black") +
geom_hline(yintercept = mod_S0$b, lty = 1, lwd = 0.75, colour = "black") +
geom_hline(yintercept = mod_S0$ci.ub, lty = 3, lwd = 0.75, colour = "black") +
geom_pointrange(aes(x = left_out, y = est, ymin = lower, ymax = upper)) +
xlab("Study left out") +
ylab("lnRR, 95% CI") +
coord_flip() +
theme(panel.grid.minor = element_blank())+
theme_bw() + theme(panel.grid.major = element_blank()) +
theme(panel.grid.minor.x = element_blank() ) +
theme(axis.text.y = element_text(size = 6))
leaveoneout_S
dat$Study_ID <- as.integer(dat$Study_ID)Interaction of stress and EE
Intercept
Enriched and stress animals are better at learning and memory.
mod_ES0 <- rma.mv(yi = lnRR_ESa, V = VCV_ES, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat)
summary(mod_ES0) ##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -50.9607 101.9214 109.9214 119.9648 110.3865
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0279 0.1669 30 no Study_ID
## sigma^2.2 0.0311 0.1765 92 no ES_ID
## sigma^2.3 0.0048 0.0690 6 no Strain
##
## Test for Heterogeneity:
## Q(df = 91) = 307.1213, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## 0.1578 0.0678 2.3273 91 0.0222 0.0231 0.2925 *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1i2_ml(mod_ES0) ## I2_total I2_Study_ID I2_ES_ID I2_Strain
## 0.82109711 0.35875892 0.40100786 0.06133033orchard_plot(mod_ES0, mod = "Int", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 24), # change font sizes
legend.title = element_text(size = 15),
legend.text = element_text(size = 13)) 
Intercept with outlier removed
dat_outliers <- dat %>%
filter(ES_ID != 88)
VCV_ES_outliers <- impute_covariance_matrix(vi = dat_outliers$lnRRV_E, cluster = dat$Study_ID, r = 0.5)
mod_ESoutlier <- rma.mv(yi = lnRR_ESa, V = VCV_ES_outliers, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat_outliers)
summary(mod_ES0)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -50.9607 101.9214 109.9214 119.9648 110.3865
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0279 0.1669 30 no Study_ID
## sigma^2.2 0.0311 0.1765 92 no ES_ID
## sigma^2.3 0.0048 0.0690 6 no Strain
##
## Test for Heterogeneity:
## Q(df = 91) = 307.1213, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## 0.1578 0.0678 2.3273 91 0.0222 0.0231 0.2925 *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1orchard_plot(mod_ESoutlier, mod = "Int", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 24), # change font sizes
legend.title = element_text(size = 15),
legend.text = element_text(size = 13)) 
Metaregression
Uni-moderator metaregression
Type of learning
The type of learning/memory response
VCV_ES1 <- impute_covariance_matrix(vi = dat1$lnRRV_ES, cluster = dat$Study_ID, r = 0.5)
mod_ES1 <- rma.mv(yi = lnRR_ESa, V = VCV_ES1, mod = ~Type_learning-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat1)
summary(mod_ES1)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -48.9339 97.8678 109.8678 124.7996 110.8921
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0344 0.1854 30 no Study_ID
## sigma^2.2 0.0246 0.1569 92 no ES_ID
## sigma^2.3 0.0040 0.0633 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 89) = 293.7418, p-val < .0001
##
## Test of Moderators (coefficients 1:3):
## F(df1 = 3, df2 = 89) = 3.5124, p-val = 0.0184
##
## Model Results:
##
## estimate se tval df pval ci.lb
## Type_learningConditioning 0.1900 0.0696 2.7284 89 0.0077 0.0516
## Type_learningHabituation 0.3107 0.1558 1.9946 89 0.0491 0.0012
## Type_learningRecognition 0.0256 0.0896 0.2852 89 0.7761 -0.1525
## ci.ub
## Type_learningConditioning 0.3284 **
## Type_learningHabituation 0.6202 *
## Type_learningRecognition 0.2036
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1r2_ml(mod_ES1) ## R2_marginal R2_coditional
## 0.07392779 0.94099939Learning_ES <- orchard_plot(mod_ES1, mod = "Type_learning", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
geom_point(aes(fill = name), size = 3, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
axis.text.x = element_text(size = 10), # change font sizes
axis.text.y = element_text(size = 10),
legend.title = element_text(size = 7),
legend.text = element_text(size = 7))
Learning_ES
Learning vs Memory
Is the assay broadly measuring learning or memory?
mod_ES2 <- rma.mv(yi = lnRR_ESa, V = VCV_ES, mod = ~Learning_vs_memory-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat)
summary(mod_ES2) ##
## Multivariate Meta-Analysis Model (k = 85; method: REML)
##
## logLik Deviance AIC BIC AICc
## -49.8593 99.7187 109.7187 121.8129 110.4979
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0282 0.1678 30 no Study_ID
## sigma^2.2 0.0338 0.1837 85 no ES_ID
## sigma^2.3 0.0052 0.0721 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 83) = 302.1628, p-val < .0001
##
## Test of Moderators (coefficients 1:2):
## F(df1 = 2, df2 = 83) = 3.0108, p-val = 0.0547
##
## Model Results:
##
## estimate se tval df pval ci.lb
## Learning_vs_memoryLearning 0.2044 0.0845 2.4175 83 0.0178 0.0362
## Learning_vs_memoryMemory 0.1379 0.0719 1.9174 83 0.0586 -0.0051
## ci.ub
## Learning_vs_memoryLearning 0.3725 *
## Learning_vs_memoryMemory 0.2810 .
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1r2_ml(mod_ES2) ## R2_marginal R2_coditional
## 0.01448953 0.92362134LvsM_ES <- orchard_plot(mod_ES2, mod = "Learning_vs_memory", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
geom_point(aes(fill = name), size = 3, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
axis.text.x = element_text(size = 10), # change font sizes
axis.text.y = element_text(size = 10),
legend.title = element_text(size = 7),
legend.text = element_text(size = 7))
LvsM_ES 
Type of reinforcement
The type of conditioning used
VCV_ES2 <- impute_covariance_matrix(vi = dat2$lnRRV_ES, cluster = dat2$Study_ID, r = 0.5)
mod_ES3 <- rma.mv(yi = lnRR_ESa, V = VCV_ES2, mod = ~ Type_reinforcement-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat2)
summary(mod_ES3)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -49.6471 99.2943 111.2943 126.2261 112.3186
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0335 0.1831 30 no Study_ID
## sigma^2.2 0.0275 0.1658 92 no ES_ID
## sigma^2.3 0.0014 0.0372 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 89) = 288.3178, p-val < .0001
##
## Test of Moderators (coefficients 1:3):
## F(df1 = 3, df2 = 89) = 3.2144, p-val = 0.0267
##
## Model Results:
##
## estimate se tval df pval ci.lb
## Type_reinforcementAppetitive 0.1482 0.1156 1.2825 89 0.2030 -0.0814
## Type_reinforcementAversive 0.1909 0.0667 2.8635 89 0.0052 0.0584
## Type_reinforcementNot applicable 0.0578 0.0798 0.7247 89 0.4705 -0.1007
## ci.ub
## Type_reinforcementAppetitive 0.3779
## Type_reinforcementAversive 0.3234 **
## Type_reinforcementNot applicable 0.2163
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1r2_ml(mod_ES3) ## R2_marginal R2_coditional
## 0.04711742 0.97883321Reinforcement_ES <- orchard_plot(mod_ES3, mod = "Type_reinforcement", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
geom_point(aes(fill = name), size = 3, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
axis.text.x = element_text(size = 10), # change font sizes
axis.text.y = element_text(size = 10),
legend.title = element_text(size = 7),
legend.text = element_text(size = 7))
Reinforcement_ES 
Type of stress
The type of stress manipulation
VCV_ES3 <- impute_covariance_matrix(vi = dat3$lnRRV_ES, cluster = dat$Study_ID, r = 0.5)
mod_ES4 <- rma.mv(yi = lnRR_ESa, V = VCV_ES3, mod = ~Type_stress_exposure-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat3)
summary(mod_ES4)##
## Multivariate Meta-Analysis Model (k = 85; method: REML)
##
## logLik Deviance AIC BIC AICc
## -45.8298 91.6596 105.6596 122.4207 107.1939
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0559 0.2365 25 no Study_ID
## sigma^2.2 0.0320 0.1788 85 no ES_ID
## sigma^2.3 0.0221 0.1487 4 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 81) = 297.6210, p-val < .0001
##
## Test of Moderators (coefficients 1:4):
## F(df1 = 4, df2 = 81) = 0.4917, p-val = 0.7418
##
## Model Results:
##
## estimate se tval df pval ci.lb
## Type_stress_exposureCombination 0.1140 0.1773 0.6429 81 0.5221 -0.2387
## Type_stress_exposureMS 0.1571 0.1392 1.1287 81 0.2624 -0.1198
## Type_stress_exposureNoise 0.2202 0.2133 1.0321 81 0.3051 -0.2043
## Type_stress_exposureRestraint 0.1788 0.1554 1.1508 81 0.2532 -0.1303
## ci.ub
## Type_stress_exposureCombination 0.4667
## Type_stress_exposureMS 0.4340
## Type_stress_exposureNoise 0.6447
## Type_stress_exposureRestraint 0.4880
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1r2_ml(mod_ES4)## R2_marginal R2_coditional
## 0.008276802 0.800570495Stressor_ES <- orchard_plot(mod_ES4, mod = "Type_stress_exposure", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
geom_point(aes(fill = name), size = 3, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
axis.text.x = element_text(size = 10), # change font sizes
axis.text.y = element_text(size = 10),
legend.title = element_text(size = 7),
legend.text = element_text(size = 7))
Stressor_ES 
Age of stress
The age at which the individuals were exposed to the stressor.
mod_ES5 <-rma.mv(yi = lnRR_ESa, V = VCV_ES, mod = ~Age_stress_exposure-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat)
summary(mod_ES5) ##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -48.8724 97.7449 111.7449 129.0862 113.1449
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0221 0.1486 30 no Study_ID
## sigma^2.2 0.0315 0.1774 92 no ES_ID
## sigma^2.3 0.0186 0.1365 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 88) = 279.7710, p-val < .0001
##
## Test of Moderators (coefficients 1:4):
## F(df1 = 4, df2 = 88) = 1.7577, p-val = 0.1446
##
## Model Results:
##
## estimate se tval df pval
## Age_stress_exposureAdolescent 0.0288 0.2073 0.1388 88 0.8899
## Age_stress_exposureAdult 0.2096 0.1331 1.5750 88 0.1188
## Age_stress_exposureEarly postnatal 0.1402 0.1098 1.2770 88 0.2050
## Age_stress_exposurePrenatal 0.3790 0.1455 2.6059 88 0.0108
## ci.lb ci.ub
## Age_stress_exposureAdolescent -0.3832 0.4408
## Age_stress_exposureAdult -0.0549 0.4741
## Age_stress_exposureEarly postnatal -0.0780 0.3583
## Age_stress_exposurePrenatal 0.0900 0.6681 *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1r2_ml(mod_ES5) ## R2_marginal R2_coditional
## 0.09585233 0.76647791Age_stress_ES<-orchard_plot(mod_ES5, mod = "Age_stress_exposure", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
geom_point(aes(fill = name), size = 3, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
axis.text.x = element_text(size = 10), # change font sizes
axis.text.y = element_text(size = 10),
legend.title = element_text(size = 7),
legend.text = element_text(size = 7))
Age_stress_ES
Stress duration
How long was the stress applied for (chronic = every day for 7 days or more)? This has the highest marginal R2 (currentl nearly 43%) - need to redo without outlier
VCV_ES4 <- impute_covariance_matrix(vi = dat4$lnRRV_ES, cluster = dat4$Study_ID, r = 0.5)
mod_ES6 <-rma.mv(yi = lnRR_ESa, V = VCV_ES4, mod = ~Stress_duration-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat4)
summary(mod_ES6) ##
## Multivariate Meta-Analysis Model (k = 89; method: REML)
##
## logLik Deviance AIC BIC AICc
## -46.2033 92.4066 102.4066 114.7362 103.1474
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0136 0.1167 29 no Study_ID
## sigma^2.2 0.0351 0.1874 89 no ES_ID
## sigma^2.3 0.0104 0.1020 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 87) = 286.8908, p-val < .0001
##
## Test of Moderators (coefficients 1:2):
## F(df1 = 2, df2 = 87) = 4.5460, p-val = 0.0132
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## Stress_durationAcute -0.0266 0.1112 -0.2388 87 0.8118 -0.2476 0.1945
## Stress_durationChronic 0.2220 0.0828 2.6815 87 0.0088 0.0575 0.3866
##
## Stress_durationAcute
## Stress_durationChronic **
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1r2_ml(mod_ES6) ## R2_marginal R2_coditional
## 0.1600392 0.8521674Duration_ES<- orchard_plot(mod_ES6, mod = "Stress_duration", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
geom_point(aes(fill = name), size = 3, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
axis.text.x = element_text(size = 10), # change font sizes
axis.text.y = element_text(size = 10),
legend.title = element_text(size = 7),
legend.text = element_text(size = 7))
Duration_ES
Social enrichment
Does EE also include a manipulation of social environment (i.e., number of individuals in EE relative to control)?
VCV_ES5 <- impute_covariance_matrix(vi = dat5$lnRRV_ES, cluster = dat5$Study_ID, r = 0.5)
mod_ES7<- rma.mv(yi = lnRR_ESa, V = VCV_ES5, mod = ~EE_social-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat5)
summary(mod_ES7)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -50.3334 100.6668 110.6668 123.1659 111.3811
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0285 0.1689 30 no Study_ID
## sigma^2.2 0.0312 0.1765 92 no ES_ID
## sigma^2.3 0.0073 0.0857 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 90) = 307.1089, p-val < .0001
##
## Test of Moderators (coefficients 1:2):
## F(df1 = 2, df2 = 90) = 2.7595, p-val = 0.0687
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## EE_socialNon-social 0.1099 0.0910 1.2080 90 0.2302 -0.0708 0.2906
## EE_socialSocial 0.2070 0.0892 2.3198 90 0.0226 0.0297 0.3843 *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1r2_ml(mod_ES7) ## R2_marginal R2_coditional
## 0.03307114 0.89402837Social_ES <- orchard_plot(mod_ES7, mod = "EE_social", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
geom_point(aes(fill = name), size = 3, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
axis.text.x = element_text(size = 10), # change font sizes
axis.text.y = element_text(size = 10),
legend.title = element_text(size = 7),
legend.text = element_text(size = 7))
Social_ES
Exercise enrichment
Does the form of enrichment include exercise through a running wheel or treadmill?
mod_ES8<- rma.mv(yi = lnRR_ESa, V = VCV_ES, mod = ~EE_exercise-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat)
summary(mod_ES8)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -50.3523 100.7046 110.7046 123.2036 111.4189
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0244 0.1561 30 no Study_ID
## sigma^2.2 0.0316 0.1778 92 no ES_ID
## sigma^2.3 0.0113 0.1062 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 90) = 294.4340, p-val < .0001
##
## Test of Moderators (coefficients 1:2):
## F(df1 = 2, df2 = 90) = 2.4502, p-val = 0.0920
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## EE_exerciseExercise 0.1298 0.0884 1.4687 90 0.1454 -0.0458 0.3053
## EE_exerciseNo exercise 0.2313 0.1085 2.1307 90 0.0358 0.0156 0.4469
##
## EE_exerciseExercise
## EE_exerciseNo exercise *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1r2_ml(mod_ES8) ## R2_marginal R2_coditional
## 0.03394203 0.83804575Exercise_ES <- orchard_plot(mod_ES8, mod = "EE_exercise", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
geom_point(aes(fill = name), size = 3, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
axis.text.x = element_text(size = 10), # change font sizes
axis.text.y = element_text(size = 10),
legend.title = element_text(size = 7),
legend.text = element_text(size = 7))
Exercise_ES
Order to treatment exposure
What order were animals exposed to stress and EE
mod_ES9 <- rma.mv(yi = lnRR_ESa, V = VCV_ES, mod = ~Exposure_order -1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat)
summary(mod_ES9)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -48.5795 97.1590 109.1590 124.0908 110.1834
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0309 0.1758 30 no Study_ID
## sigma^2.2 0.0303 0.1741 92 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 89) = 295.5040, p-val < .0001
##
## Test of Moderators (coefficients 1:3):
## F(df1 = 3, df2 = 89) = 4.1173, p-val = 0.0088
##
## Model Results:
##
## estimate se tval df pval ci.lb
## Exposure_orderConcurrently 0.2255 0.1360 1.6583 89 0.1008 -0.0447
## Exposure_orderEnrichment first -0.2430 0.2047 -1.1871 89 0.2384 -0.6496
## Exposure_orderStress first 0.1597 0.0558 2.8623 89 0.0052 0.0488
## ci.ub
## Exposure_orderConcurrently 0.4958
## Exposure_orderEnrichment first 0.1637
## Exposure_orderStress first 0.2705 **
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1r2_ml(mod_ES9)## R2_marginal R2_coditional
## 0.1525456 1.0000000Order_ES <- orchard_plot(mod_ES9, mod = "Exposure_order", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
geom_point(aes(fill = name), size = 3, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
axis.text.x = element_text(size = 10), # change font sizes
axis.text.y = element_text(size = 10),
legend.title = element_text(size = 7),
legend.text = element_text(size = 7))
Order_ES 
Age of enrichment
What age were individuals exposed to EE
VCV_ES6 <- impute_covariance_matrix(vi = dat6$lnRRV_ES, cluster = dat6$Study_ID, r = 0.5)
mod_ES10 <- rma.mv(yi = lnRR_ESa, V = VCV_ES6, mod = ~Age_EE_exposure-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat6)
summary(mod_ES10) ##
## Multivariate Meta-Analysis Model (k = 88; method: REML)
##
## logLik Deviance AIC BIC AICc
## -45.0007 90.0014 100.0014 112.2731 100.7514
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0315 0.1776 29 no Study_ID
## sigma^2.2 0.0277 0.1663 88 no ES_ID
## sigma^2.3 0.0022 0.0467 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 86) = 291.6264, p-val < .0001
##
## Test of Moderators (coefficients 1:2):
## F(df1 = 2, df2 = 86) = 3.5102, p-val = 0.0342
##
## Model Results:
##
## estimate se tval df pval ci.lb
## Age_EE_exposureAdolescent 0.1616 0.0666 2.4263 86 0.0173 0.0292
## Age_EE_exposureAdult 0.1527 0.1041 1.4661 86 0.1463 -0.0543
## ci.ub
## Age_EE_exposureAdolescent 0.2941 *
## Age_EE_exposureAdult 0.3597
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1r2_ml(mod_ES10) ## R2_marginal R2_coditional
## 0.0002073979 0.9645265244Age_enrichment_ES <- orchard_plot(mod_ES10, mod = "Age_EE_exposure", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
geom_point(aes(fill = name), size = 3, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
axis.text.x = element_text(size = 10), # change font sizes
axis.text.y = element_text(size = 10),
legend.title = element_text(size = 7),
legend.text = element_text(size = 7))
Age_enrichment_ES
Multi-moderator fullmodel
dat_ESfm <- dat %>%
filter(Type_learning %in% c("Recognition", "Habituation", "Conditioning"),
Type_reinforcement %in% c("Appetitive", "Aversive", "Not applicable"),
EE_social %in% c("Social", "Non-social"),
Age_EE_exposure %in% c("Adult", "Adolescent"),
Type_stress_exposure %in% c("Restraint", "Noise", "MS", "Combination"),
Stress_duration %in% c("Chronic", "Acute"),
Age_stress_exposure %in% c("Prenatal", "Early postnatal", "Adult"))
VCV_ESfm <- impute_covariance_matrix(vi = dat_ESfm$lnRRV_ES, cluster = dat$Study_ID, r = 0.5)
mod_ESfm <- rma.mv(yi = lnRR_Sa, V = VCV_ESfm, mod = ~Type_learning-1 + Learning_vs_memory-1 + Type_reinforcement-1 + EE_social-1 + EE_exercise-1 + Age_EE_exposure-1 + Type_stress_exposure-1 + Age_stress_exposure-1 + Stress_duration-1 + Exposure_order, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat_ESfm)
#summary(mod_ESfm)
#r2_ml(mod_ESfm)
res_ESfm <- dredge(mod_ESfm, trace=2)##
|
| | 0%
|
| | 1%
|
|= | 1%
|
|= | 2%
|
|== | 2%
|
|== | 3%
|
|== | 4%
|
|=== | 4%
|
|=== | 5%
|
|==== | 5%
|
|==== | 6%
|
|===== | 6%
|
|===== | 7%
|
|===== | 8%
|
|====== | 8%
|
|====== | 9%
|
|======= | 9%
|
|======= | 10%
|
|======= | 11%
|
|======== | 11%
|
|======== | 12%
|
|========= | 12%
|
|========= | 13%
|
|========== | 14%
|
|========== | 15%
|
|=========== | 15%
|
|=========== | 16%
|
|============ | 17%
|
|============ | 18%
|
|============= | 18%
|
|============= | 19%
|
|============== | 19%
|
|============== | 20%
|
|============== | 21%
|
|=============== | 21%
|
|=============== | 22%
|
|================ | 22%
|
|================ | 23%
|
|================ | 24%
|
|================= | 24%
|
|================= | 25%
|
|================== | 25%
|
|================== | 26%
|
|=================== | 26%
|
|=================== | 27%
|
|=================== | 28%
|
|==================== | 28%
|
|==================== | 29%
|
|===================== | 29%
|
|===================== | 30%
|
|===================== | 31%
|
|====================== | 31%
|
|====================== | 32%
|
|======================= | 32%
|
|======================= | 33%
|
|======================== | 34%
|
|======================== | 35%
|
|========================= | 35%
|
|========================= | 36%
|
|========================== | 37%
|
|========================== | 38%
|
|=========================== | 38%
|
|=========================== | 39%
|
|============================ | 39%
|
|============================ | 40%
|
|============================ | 41%
|
|============================= | 41%
|
|============================= | 42%
|
|============================== | 42%
|
|============================== | 43%
|
|============================== | 44%
|
|=============================== | 44%
|
|=============================== | 45%
|
|================================ | 45%
|
|================================ | 46%
|
|================================= | 46%
|
|================================= | 47%
|
|================================= | 48%
|
|================================== | 48%
|
|================================== | 49%
|
|=================================== | 49%
|
|=================================== | 50%
|
|=================================== | 51%
|
|==================================== | 51%
|
|==================================== | 52%
|
|===================================== | 52%
|
|===================================== | 53%
|
|===================================== | 54%
|
|====================================== | 54%
|
|====================================== | 55%
|
|======================================= | 55%
|
|======================================= | 56%
|
|======================================== | 56%
|
|======================================== | 57%
|
|======================================== | 58%
|
|========================================= | 58%
|
|========================================= | 59%
|
|========================================== | 59%
|
|========================================== | 60%
|
|========================================== | 61%
|
|=========================================== | 61%
|
|=========================================== | 62%
|
|============================================ | 62%
|
|============================================ | 63%
|
|============================================= | 64%
|
|============================================= | 65%
|
|============================================== | 65%
|
|============================================== | 66%
|
|=============================================== | 67%
|
|=============================================== | 68%
|
|================================================ | 68%
|
|================================================ | 69%
|
|================================================= | 69%
|
|================================================= | 70%
|
|================================================= | 71%
|
|================================================== | 71%
|
|================================================== | 72%
|
|=================================================== | 72%
|
|=================================================== | 73%
|
|=================================================== | 74%
|
|==================================================== | 74%
|
|==================================================== | 75%
|
|===================================================== | 75%
|
|===================================================== | 76%
|
|====================================================== | 76%
|
|====================================================== | 77%
|
|====================================================== | 78%
|
|======================================================= | 78%
|
|======================================================= | 79%
|
|======================================================== | 79%
|
|======================================================== | 80%
|
|======================================================== | 81%
|
|========================================================= | 81%
|
|========================================================= | 82%
|
|========================================================== | 82%
|
|========================================================== | 83%
|
|=========================================================== | 84%
|
|=========================================================== | 85%
|
|============================================================ | 85%
|
|============================================================ | 86%
|
|============================================================= | 87%
|
|============================================================= | 88%
|
|============================================================== | 88%
|
|============================================================== | 89%
|
|=============================================================== | 89%
|
|=============================================================== | 90%
|
|=============================================================== | 91%
|
|================================================================ | 91%
|
|================================================================ | 92%
|
|================================================================= | 92%
|
|================================================================= | 93%
|
|================================================================= | 94%
|
|================================================================== | 94%
|
|================================================================== | 95%
|
|=================================================================== | 95%
|
|=================================================================== | 96%
|
|==================================================================== | 96%
|
|==================================================================== | 97%
|
|==================================================================== | 98%
|
|===================================================================== | 98%
|
|===================================================================== | 99%
|
|======================================================================| 99%
|
|======================================================================| 100%res_ESfm2<- subset(res_ESfm, delta <= 6, recalc.weights=FALSE)
importance(res_ESfm2) ## Learning_vs_memory Stress_duration Age_EE_exposure
## Sum of weights: 0.81 0.76 0.33
## N containing models: 44 37 20
## EE_social Age_stress_exposure EE_exercise
## Sum of weights: 0.24 0.20 0.19
## N containing models: 17 10 13
## Type_stress_exposure Type_learning Exposure_order
## Sum of weights: 0.18 0.10 0.05
## N containing models: 12 10 4
## Type_reinforcement
## Sum of weights: 0.01
## N containing models: 2Publication bias & sensitivity analysis
Funnel_ES<-funnel(mod_ESfm, xlab = "lnRR", ylab = "Standard Error")
#year published was scaled previously under stress PB
dat_ESfm$sqrt_inv_e_n <- with(dat_ESfm, sqrt(1/CC_n + 1/EC_n + 1/ES_n + 1/CS_n))
PB_MR_ES<- rma.mv(lnRR_ESa, lnRRV_ES, mods = ~1 + sqrt_inv_e_n + Year_published + Type_learning + Type_reinforcement + EE_social + EE_exercise + Age_stress_exposure, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain), method = "REML", test = "t",
data = dat_ESfm)
estimates_PB_MR_ES<- estimates.CI(PB_MR_ES)
estimates_PB_MR_ES| estimate | mean | lower | upper |
|---|---|---|---|
| intrcpt | 44.1112313 | -25.6589074 | 113.8813700 |
| sqrt_inv_e_n | -0.6897056 | -1.8750276 | 0.4956163 |
| Year_published | -0.0216674 | -0.0562695 | 0.0129347 |
| Type_learningHabituation | 0.2280781 | -0.4019047 | 0.8580609 |
| Type_learningRecognition | -0.0104748 | -0.5598110 | 0.5388614 |
| Type_reinforcementAversive | 0.0003694 | -0.4348552 | 0.4355939 |
| Type_reinforcementNot applicable | -0.2293071 | -0.9140504 | 0.4554362 |
| EE_socialSocial | 0.1166764 | -0.2532547 | 0.4866075 |
| EE_exerciseNo exercise | 0.1207767 | -0.2456524 | 0.4872057 |
| Age_stress_exposureEarly postnatal | 0.1583410 | -0.2028311 | 0.5195131 |
| Age_stress_exposurePrenatal | 0.1954658 | -0.3134537 | 0.7043854 |
#no effect of inv_effective sample size or year published
#leave-one-out analysis
dat$Study_ID <- as.factor(dat$Study_ID)
LeaveOneOut_effectsize <- list()
for(i in 1:length(levels(dat$Study_ID))){
LeaveOneOut_effectsize[[i]] <- rma.mv(yi = lnRR_Ea, V = lnRRV_E,
random = list(~1 | Study_ID,~1| ES_ID, ~1 | Strain),
method = "REML", data = dat[dat$Study_ID
!= levels(dat$Study_ID)[i], ])}
# writing function for extracting est, ci.lb, and ci.ub from all models
est.func <- function(mod_E0){
df <- data.frame(est = mod_E0$b, lower = mod_E0$ci.lb, upper = mod_E0$ci.ub)
return(df)
}
#using dplyr to form data frame
MA_CVR <- lapply(LeaveOneOut_effectsize, function(x) est.func(x))%>% bind_rows %>% mutate(left_out = levels(dat$Study_ID))
#telling ggplot to stop reordering factors
MA_CVR$left_out<- as.factor(MA_CVR$left_out)
MA_CVR$left_out<-factor(MA_CVR$left_out, levels = MA_CVR$left_out)
#plotting
leaveoneout_ES <- ggplot(MA_CVR) +
geom_hline(yintercept = 0, lty = 2, lwd = 1) +
geom_hline(yintercept = mod_E0$ci.lb, lty = 3, lwd = 0.75, colour = "black") +
geom_hline(yintercept = mod_E0$b, lty = 1, lwd = 0.75, colour = "black") +
geom_hline(yintercept = mod_E0$ci.ub, lty = 3, lwd = 0.75, colour = "black") +
geom_pointrange(aes(x = left_out, y = est, ymin = lower, ymax = upper)) +
xlab("Study left out") +
ylab("lnRR, 95% CI") +
coord_flip() +
theme(panel.grid.minor = element_blank())+
theme_bw() + theme(panel.grid.major = element_blank()) +
theme(panel.grid.minor.x = element_blank() ) +
theme(axis.text.y = element_text(size = 6))
leaveoneout_ES
dat$Study_ID <- as.integer(dat$Study_ID)Combined orchard plot
mod_list1 <- list(mod_E0, mod_S0, mod_ES0)
mod_res1 <- lapply(mod_list1, function(x) mod_results(x, mod = "Int"))
merged1 <- submerge(mod_res1[[3]], mod_res1[[2]], mod_res1[[1]], mix = T)
merged1$mod_table$name <- factor(merged1$mod_table$name, levels = c("Intrcpt1",
"Intrcpt2", "Intrcpt3"),
labels = rev(c("Enrichment ME", "Stress ME", "Interaction")))
merged1$data$moderator <- factor(merged1$data$moderator, levels = c("Intrcpt1",
"Intrcpt2", "Intrcpt3"),
labels = rev(c("Enrichment ME", "Stress ME", "Interaction")))
orchard1<- orchard_plot(merged1, mod = "Int", xlab = "lnRR", angle = 0) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
xlim(-2,4.5) +
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 15), # change font sizes
legend.title = element_text(size = 10),
legend.text = element_text(size = 10))
orchard1
Age of enrichment
The age at which the individuals were exposed to environmental enrichment.
This needs to be redone after recategorising
mod_ES9 <- rma.mv(yi = lnRR_ESa, V = lnRRV_ES, mod = ~Age_EE_exposure-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat)
summary(mod_ES9) ##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -47.5384 95.0767 107.0767 122.0086 108.1011
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0450 0.2122 30 no Study_ID
## sigma^2.2 0.0194 0.1393 92 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 89) = 302.2511, p-val < .0001
##
## Test of Moderators (coefficients 1:3):
## F(df1 = 3, df2 = 89) = 4.6525, p-val = 0.0046
##
## Model Results:
##
## estimate se tval df pval ci.lb
## Age_EE_exposureAdolescent 0.2006 0.0599 3.3519 89 0.0012 0.0817
## Age_EE_exposureAdult 0.1526 0.0979 1.5582 89 0.1227 -0.0420
## Age_EE_exposureEarly postnatal -0.1674 0.3086 -0.5424 89 0.5889 -0.7805
## ci.ub
## Age_EE_exposureAdolescent 0.3196 **
## Age_EE_exposureAdult 0.3472
## Age_EE_exposureEarly postnatal 0.4457
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1r2_ml(mod_ES9) ## R2_marginal R2_coditional
## 0.082042 1.000000# Orchard plot
orchard_plot(mod_ES9, mod = "Age_EE_exposure", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 24), # change font sizes
legend.title = element_text(size = 15),
legend.text = element_text(size = 13)) 
‘Pairwise’ effect sizes
Enrichment relative to control
VCV_E20 <- impute_covariance_matrix(vi = dat$lnRRV_E2, cluster = dat$Study_ID, r = 0.5)
#Model doesn't converge with VCV
mod_E20 <- rma.mv(yi = lnRR_E2a, V = VCV_E20, random = list(~1|Study_ID,
~ 1|Strain,
~1|ES_ID),
test = "t",
data = dat,
control=list(optimizer="optim", optmethod="Nelder-Mead"))
summary(mod_E20) ##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -80.4236 160.8471 168.8471 178.8905 169.3122
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0209 0.1444 30 no Study_ID
## sigma^2.2 0.0000 0.0001 6 no Strain
## sigma^2.3 0.0000 0.0000 92 no ES_ID
##
## Test for Heterogeneity:
## Q(df = 91) = 23.4395, p-val = 1.0000
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## 0.1376 0.0729 1.8860 91 0.0625 -0.0073 0.2825 .
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1i2_ml(mod_E20) ## I2_total I2_Study_ID I2_Strain I2_ES_ID
## 8.461245e-02 8.461241e-02 4.008740e-08 0.000000e+00orchard_plot(mod_E20, mod = "Int", xlab = "lnRR")
Stress relative to control
VCV_S20 <- impute_covariance_matrix(vi = dat$lnRRV_S2, cluster = dat$Study_ID, r = 0.5)
mod_S20 <- rma.mv(yi = lnRR_S2a, V = VCV_S20, random = list(~1|Study_ID,
~ 1|Strain,
~1|ES_ID),
test = "t",
data = dat)
summary(mod_S20) ##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -89.4392 178.8785 186.8785 196.9219 187.3436
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0300 0.1732 30 no Study_ID
## sigma^2.2 0.0000 0.0000 6 no Strain
## sigma^2.3 0.0000 0.0000 92 no ES_ID
##
## Test for Heterogeneity:
## Q(df = 91) = 37.2472, p-val = 1.0000
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## -0.0735 0.0790 -0.9307 91 0.3545 -0.2304 0.0834
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1i2_ml(mod_S20) ## I2_total I2_Study_ID I2_Strain I2_ES_ID
## 1.143516e-01 1.143515e-01 7.749001e-09 0.000000e+00orchard_plot(mod_S20, mod = "Int", xlab = "lnRR")
Enrichment + stress relative to control
VCV_ES20 <- impute_covariance_matrix(vi = dat$lnRRV_ES2, cluster = dat$Study_ID, r = 0.5)
mod_ES20 <- rma.mv(yi = lnRR_ES2a, V = VCV_ES20, random = list(~1|Study_ID,
~ 1|Strain,
~1|ES_ID),
test = "t",
data = dat)
summary(mod_ES20) ##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -86.3233 172.6465 180.6465 190.6900 181.1116
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0238 0.1541 30 no Study_ID
## sigma^2.2 0.0000 0.0000 6 no Strain
## sigma^2.3 0.0000 0.0000 92 no ES_ID
##
## Test for Heterogeneity:
## Q(df = 91) = 29.6180, p-val = 1.0000
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## 0.1100 0.0772 1.4247 91 0.1577 -0.0434 0.2634
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1i2_ml(mod_ES20) ## I2_total I2_Study_ID I2_Strain I2_ES_ID
## 8.967119e-02 8.967119e-02 4.077839e-10 1.167794e-10orchard_plot(mod_ES20, mod = "Int", xlab = "lnRR")
Enrichment + stress relative to stress
VCV_E30 <- impute_covariance_matrix(vi = dat$lnRRV_E3, cluster = dat$Study_ID, r = 0.5)
mod_E30 <- rma.mv(yi = lnRR_E3a, V = VCV_E30, random = list(~1|Study_ID,
~ 1|Strain,
~1|ES_ID),
test = "t",
data = dat)
summary(mod_E30) ##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -90.8497 181.6995 189.6995 199.7429 190.1646
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0097 0.0986 30 no Study_ID
## sigma^2.2 0.0000 0.0001 6 no Strain
## sigma^2.3 0.0000 0.0000 92 no ES_ID
##
## Test for Heterogeneity:
## Q(df = 91) = 27.4127, p-val = 1.0000
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## 0.1293 0.0674 1.9192 91 0.0581 -0.0045 0.2632 .
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1i2_ml(mod_E30) ## I2_total I2_Study_ID I2_Strain I2_ES_ID
## 4.001741e-02 4.001738e-02 2.945906e-08 1.186013e-14orchard_plot(mod_E30, mod = "Int", xlab = "lnRR")
Enrichment + stress relative to enrichment
VCV_S30 <- impute_covariance_matrix(vi = dat$lnRRV_S3, cluster = dat$Study_ID, r = 0.5)
mod_S30 <- rma.mv(yi = lnRR_S3a, V = VCV_S30, random = list(~1|Study_ID,
~ 1|Strain,
~1|ES_ID),
test = "t",
data = dat,
control=list(optimizer="optim", optmethod="Nelder-Mead"))
summary(mod_S30) ##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -80.0114 160.0229 168.0229 178.0663 168.4880
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0000 0.0000 30 no Study_ID
## sigma^2.2 0.0000 0.0000 6 no Strain
## sigma^2.3 0.0000 0.0000 92 no ES_ID
##
## Test for Heterogeneity:
## Q(df = 91) = 15.4501, p-val = 1.0000
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## -0.0120 0.0618 -0.1942 91 0.8465 -0.1347 0.1107
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1i2_ml(mod_S30) ## I2_total I2_Study_ID I2_Strain I2_ES_ID
## 0 0 0 0orchard_plot(mod_S30, mod = "Int", xlab = "lnRR")
Combined orchard plot
mod_list2 <- list(mod_S30, mod_E30, mod_ES20, mod_S20, mod_E20) #rearranged the order so that it matches intext results
mod_res2 <- lapply(mod_list2, function(x) mod_results(x, mod = "Int"))
merged2 <- submerge(mod_res2[[1]], mod_res2[[2]], mod_res2[[3]], mod_res2[[4]], mod_res2[[5]], mix = T)
merged2$mod_table$name <- factor(merged2$mod_table$name, levels = c("Intrcpt1",
"Intrcpt2", "Intrcpt3", "Intrcpt4", "Intrcpt5"),
labels = rev(c("EC/CC", "CS/CC", "ES/CC", "ES/CS", "ES/EC")))
merged2$data$moderator <- factor(merged2$data$moderator, levels = c("Intrcpt1",
"Intrcpt2", "Intrcpt3", "Intrcpt4", "Intrcpt5"),
labels = rev(c("EC/CC", "CS/CC", "ES/CC", "ES/CS", "ES/EC")))
orchard2 <- orchard_plot(merged2, mod = "Int", xlab = "lnRR", angle = 0) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
xlim(-2,4.5) +
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 15), # change font sizes
legend.title = element_text(size = 10),
legend.text = element_text(size = 10))
orchard2
Panel of ‘focal’ ES and ‘pairwise’ ES orchard plots
p1 <- orchard1 + orchard2 + plot_annotation(tag_levels = 'A')
p1
Panel of metaregressions
Environmental enrichment
#Enrichment
E_mod <- (LvsM_E + Learning_E + Reinforcement_E)/ (Age_E + Exercise_E + Social_E) + plot_annotation(tag_levels = 'A')
E_mod
Stress
S_mod <- (LvsM_S + Learning_S + Reinforcement_S) / (Age_S + Stressor + Duration_S) + plot_annotation(tag_levels = 'A')
S_mod
Interaction
ES_mod <- plot_grid(LvsM_ES, Learning_ES, Reinforcement_ES, Age_enrichment_ES, Age_stress_ES, Order_ES, Exercise_ES, Social_ES, Stressor_ES, Duration_ES,
labels = "AUTO", ncol = 5)
ES_mod
Modelling with SMD
Environmental Enrichment
Intercept model
Why does the funnel plot look so strange and why is I2 higher?
mod_E0a <- rma.mv(yi = SMD_Ea, V = VCV_Ea, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat)
summary(mod_E0a)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -113.1872 226.3745 234.3745 244.4179 234.8396
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.3274 0.5722 30 no Study_ID
## sigma^2.2 0.4890 0.6993 92 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Heterogeneity:
## Q(df = 91) = 12552.8093, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## 0.7290 0.1333 5.4689 91 <.0001 0.4642 0.9938 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1i2_ml(mod_E0a)## I2_total I2_Study_ID I2_ES_ID I2_Strain
## 9.972502e-01 3.999606e-01 5.972896e-01 1.643432e-09Stress
Intercept model
Why does the funnel plot look so strange and why is I2 higher?
mod_S0a <- rma.mv(yi = SMD_Sa, V = VCV_Sa, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat)
summary(mod_S0a) ##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -121.4344 242.8688 250.8688 260.9122 251.3339
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.5602 0.7484 30 no Study_ID
## sigma^2.2 0.5409 0.7355 92 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Heterogeneity:
## Q(df = 91) = 28682.7938, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## -0.4360 0.1627 -2.6795 91 0.0088 -0.7592 -0.1128 **
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1i2_ml(mod_S0a) ## I2_total I2_Study_ID I2_ES_ID I2_Strain
## 9.979596e-01 5.076959e-01 4.902636e-01 2.171066e-09Interaction
Intercept model
Why does the funnel plot look so strange and why is I2 higher?
mod_ES0a <- rma.mv(yi = SMD_ESa, V = VCV_ESa, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat)
summary(mod_ES0a)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -126.8085 253.6170 261.6170 271.6604 262.0821
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.7153 0.8457 30 no Study_ID
## sigma^2.2 0.5861 0.7656 92 no ES_ID
## sigma^2.3 0.0000 0.0001 6 no Strain
##
## Test for Heterogeneity:
## Q(df = 91) = 11648.6783, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## 0.6961 0.1810 3.8449 91 0.0002 0.3365 1.0557 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1i2_ml(mod_ES0a) ## I2_total I2_Study_ID I2_ES_ID I2_Strain
## 9.931281e-01 5.458585e-01 4.472695e-01 1.577079e-08
Social enrichment
Does EE also include a manipulation of social environment? Note that we excluded any studies that exclusively used social enrichment.s